###Create FELM ME Function
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.05, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### defining function to get marginal effects for specific levels of the treatment variable
  getmfx_high_low <- function(betas, data, treatment, moderator, treatment_val) {
    betas[treatment] * treatment_val + betas[paste0(treatment, ":", moderator)] * data[, moderator] * treatment_val
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting  marginal effects for high and low cases of treatment variable
  data[, "marginal_effects_treatment_low"] <- getmfx_high_low(coef(regression_model), data, treatment, moderator, quantile(data[,treatment], 0.05))
  data[, "marginal_effects_treatment_high"] <- getmfx_high_low(coef(regression_model), data, treatment, moderator, quantile(data[,treatment], 0.95))
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  # Marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_hline(yintercept = 0, color = "grey70", size = 0.5) +  # light reference line at y=0
    geom_point(aes(y = marginal_effects)) + 
    geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.2) +  
    theme_fivethirtyeight() +
    theme(
      plot.title = element_text(size = 11.5, hjust = 0.5),
      axis.title = element_text(size = 10),
      panel.background = element_rect(fill = "white", color = NA),
      plot.background  = element_rect(fill = "white", color = NA),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      axis.line = element_line(color = "black")  # keep dark axis lines
    ) +
    scale_x_continuous(breaks = c(0, 1), labels = c("pre-insurgency onset", "post-insurgency onset")) +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of", treatment_translation, "\non", dependent_variable_translation))
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}